Converging toward a practical solution of the Holstein 

molecular crystal model 

A. H. Romero 1 ' 3 , David W. Brown 2 and Katja Lindenberg 3 

1 Department of Physics, 
University of California, San Diego, La Jolla, CA 92093-0354 

2 Institute for Nonlinear Science, 
University of California, San Diego, La Jolla, CA 92093-0402 

3 Department of Chemistry and Biochemistry, 
University of California, San Diego, La Jolla, CA 92093-0340 

February 1, 2008 

Abstract 

We present selected results for the Holstein molecular crystal model in one space 
dimension as determined by the Global-Local variational method, including complete 
polaron energy bands, ground state energies, and effective masses. We juxtapose our 
results with specific comparable results of numerous other methodologies of current 
interest, including quantum Monte Carlo, cluster diagonalization, dynamical mean field 
theory, density matrix renormalization group, semiclassical analysis, weak-coupling 
perturbation theory, and strong-coupling perturbation theory. Taken as a whole, these 
methodologies are mutually confirming and provide a comprehensive and quantitatively 
accurate description of polaron properties in essentially any regime. In particular, this 
comparison confirms the Global-Local variational method as being highly accurate 
over a wide range of the polaron parameter space, from the non-adiabatic limit to the 
extremes of high adiabaticity, from weak coupling through intermediate coupling to 
strong coupling. 
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1 Introduction 



Great progress has been made in the last few years toward achieving the practical solution 
of the Holstein molecular crystal model,Si a practical solution being one which yields a 
sufficiently complete spectrum of information to a sufficiently high accuracy to satisfy most 
purposes. A number of different methodologies, each important for the particular strength 
they lend to the problem, are proving increasingly mutually consistent if not always straight- 
forwardly convergent in their conclusions. In this paper, we compare a suite of our own 
recent result^ obtained using the Global-Local variational method with comparable re- 
sults of other contemporary approaches, including quantum Monte Carlo (QMC),oEj cluster 
diagonalization,E^rEZI dynamical mean field theory (DMFT),c§t2l density matrix renormaliza- 
tion group (DMRG)JMta as well as both weak-coup line perturbation theory (WCPT)S>§ 
and strong-coupling perturbation theory (SCPT).00e1 Every such comparison we have 
been able to implement has shown the Global-Local results to be among the best currently 
available, and where a better result obtains by another method, the quantitative discrepancy 
involved tends to be satisfyingly small. Thus, although it is not the case that our methods 
are necessarily the most accurate possible in every instance, our methods are valid over large 
enough a region of the polaron parameter space, and accurate enough over their region of 
validity to be sufficient to a wealth of practical purposes. 

In several recent works (10 - U@ - Illi) we have presented a sequence of increasingly 
refined variational approaches to the problem of determining the lowest polaron energy band 
and the associated energy-momentum eigenfunctions for the Holstein molecular crystal model 
in one space dimension. The Global-Local method generalizes those of ToyozawalzhEl and 
MerrifieldSli by including local electron-phonon correlation channels under-represented in 
the former and global electron-phonon correlation channels under-represented in the latter, 
both of which are of particular significance in the local structure of the polaron. While in the 
present work we lean heavily on the results of Paper III because of its superior accuracy, this 
sequential approach has proven instrumental in motivating our present study and presaging 
some of our conclusions. 

This paper is organized as follows: In Section 2, we present the model and states upon 
which the present work is based, and set down notation. In Section 3, we focus on the global 
ground state energy, displaying specific results according to our own method and comparing 
our results with those of other authors and certain approximate formulas. In Section 4, we 
present some particular examples of complete polaron energy bands and compare Global- 
Local results with specific results obtained by the DMRG method and by direct cluster 
diagonalization. In Section 5, we turn to the polaron effective mass, computing effective 
mass curves cutting swaths through the polaron parameter space in several regimes, and make 
specific comparisons of Global-Local results with those of a variety of competing approaches. 
Conclusions are summarized in Section 6. 

2 Model, States, Method 

As our system Hamiltonian, we choose the traditional Holstein Hamiltonian.01 

H = H el + H ph + H el - ph , (1) 



H el = EJ2 a l a n ~ J a i( a n+l + On-l) , 



(2) 



n n 



(3) 



n 



H el - ph = -ghu;Y,ala n (bi + bn) , 



(4) 



in which a^ n creates an electron in the rigid-lattice Wannier state at site n, and b^ n creates 
a quantum of vibrational energy in the Einstein oscillator at site n. We presume periodic 
boundary conditions on a one-dimensional lattice of N sites. The electron transfer inte- 
gral between nearest-neighbor sites is denoted by J, u is the Einstein frequency, and g is 
the dimensionless local coupling strength. (Except where displayed for formulaic clarity, 
the reference energy E is set to zero throughout.) The lattice constant does not appear 
explicitly in this formulation provided wave vectors are measured relative to it, as will be 
our convention. Two dimensionless control parameters can be constructed from the three 
principal Hamiltonian parameters in different ways. One such non-dimensionalizing scheme 
involves selecting the phonon quantum %u as the unit of energy; in these terms, the natu- 
ral dimensionless parameters are the remaining coefficients of the electronic hopping term 
(J/huj) and electron-phonon interaction term (g). This scheme is particularly appropriate 
when considering dependences on J and/or g at fixed lu, such as we shall be concerned with 
in most of this paper. This scheme is not so convenient near the adiabatic limit, however, 
where both J/hu and g diverge in a certain fixed relationship. In the adiabatic regime, it is 
convenient to non-dimensionalize by selecting the electron half-bandwidth 2J as the unit of 
energy; in these terms, the natural dimensionless parameters are the remaining coefficients 
of the phonon energy (7 = hu/2J) and the electron-phonon interaction term, expressed in 



the form y A/7, where A = g 2 Tiu/2J. (Here, we follow the convention of including the site 

coordination number z in the definition of A, such that A = g 2 hu/zJ.) The adiabatic limit 
is reached by allowing 7 to vanish at arbitrarily fixed A. We distinguish these two options 
as non-adiabatic and adiabatic scaling conventions, respectively. Except where explicitly 
noted, we conform to the non-adiabatic scaling convention. 

Our central interest in this paper is in the polaron energy band, computed as 



These states are eigenf unctions of the appropriate total momentum operator and orthogonal 
for distinct k, making variations for distinct k independent.il Since selected k values such as 
the global polaron ground state (k = 0) are generally insufficient to convey a very complete 
picture of polaron structure, we have found it important to examine complete variational 
solutions both for quality assurance and to extract the best description of polaron structure. 
By "complete" variational solutions, we mean a set of N variational energies E(k) and 





(5) 
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N polaron Bloch states the latter being described by a distinct set of variational 

parameters for each k. The set of E{k) so produced constitute an estimate (upper bound) 
for the polaron energy band.Si! 

The Global-Local method@ represents polaron structure through three classes of varia- 
tional parameters {a^,/?g, 7g}- 

I*> = E^X a -n4 a exp{-A^B^ (7) 

nn a { q ) 

We note here for later emphasis that a solution of the Global-Local method for any particular 
k is contained in the 3iV complex quantities {a%, Pg , 7^} , and that a complete solution 
is contained in 3N 2 complex quantities, reducible to 0{|iV 2 } independent real quantities 
utilizing Hamiltonian symmetries. We typically computed complete energy band structures 
on 32-site lattices for every parameter set (J/tkj, g) we considered; thus, the solution for any 
particular k is encoded in 96 independent real numbers (48 for the ground state), and the 
complete solution (all k's) in 1536 independent real numbers (see, e.g., Figure |3] in Section 
4). 

Since our solutions are encoded in a relatively small number of variational parameters 
(compared to some other calculations we shall consider), even "complete" band structure 
solutions can be stored compactly, allowing a "library" of polaron band structures to be 
archived and revisited at leisure. 

We solve the variational equations by relaxation techniques!' 0H through which an ini- 
tializing state is iteratively refined toward the self-consistent target state. Unlike many other 
methods, nearby solutions can be used to initialize new calculations, accelerating convergence 
and reducing the need to obtain new solutions from scratch. Such initializing solutions might 
be obtained from nearby points in parameter space, for example, or might be lower-precision 
solutions obtained from prior calculations at the same parameters and k. A library of po- 
laron structures sampling the polaron parameter space thus facilitates the acquisition of new 
solutions. 

The computational time required for our calculation varies according to the region of the 
parameter space examined, the error tolerance required, and the "scope" of the calculation. 
In the intermediate coupling regime, complete band structures computed to tolerances ad- 
equate for this paper take one to two hours on a single-processor Sun Microsystems Ultra 
Sparc I workstation; optimized to compute only the effective mass, only a few minutes per 
effective mass value are required. 

More intensive calculation is required as one moves out from the intermediate coupling 
region, but in no reasonably general case for which we have achieved converged results has 
a complete set of N band energies and N Bloch states required more than two to four hours 
for one (J/hu,g) point. Some special cases presented greater difficulty; e.g., the limit of 
weak coupling for J/tkv ~ 1/4 proved inaccessible due to deteriorating numerical precision, 
and convergence in the vicinity of the self-trapping transition grew more challenging with 
increasing adiabaticity. 

The self-consistency equations that follow from applying the variational principle to this 
class of trial states, the method of solving those equations, and sample results have been 
detailed in Paper III.! In the following sections, we compare polaron energy band character- 
istics as determined by the Global-Local method with a variety of alternative descriptions, 
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but apart from the one complete solution illustrated in Figure [| we defer to a subsequent 
workQ any detailed discussion of the particulars of internal polaron structure as reflected in 
the variational parameters themselves. 

3 Ground State Energy 

Without question, the single most important state in any quantum system is the ground 
state, and by association, the single most important energy is the ground state energy. 
That being said, we emphasize quickly that a particular numerical value of the ground state 
energy in itself contains very little information about the system that is of practical use. For 
most purposes, we require relationships between energies; e.g., effective masses, bandwidths, 
densities of states. For such purposes we require multiple energies, all but one of which are 
excited states in the global sense. Variational approaches such as ours divide the total space 
of the problem into subspaces of the total crystal momentum within which the ground state 
of each k sector is determined independently. In this sense, the global ground state energy 
E(0) is only one of iV independent K-sector ground state energies E{k) computed on an equal 
footing. There is no guarantee that the numerical value of E(0) constitutes any more (or 
less) accurate an estimate of its particular target value than any other E(k) so determined. 

Nonetheless, it is the global ground state energy E(0) that is the most common denom- 
inator of diverse variational approaches, including those based on localized states or other 
states that may not be well adapted to the properties expected of energy-momentum eigen- 
functions. Thus, we have compiled in Tables |I| and |2] a number of ground state energies from 
our own calculations together with a sampling of others for two particular sets of parameters 
of the Holstein Hamiltonian. 

We have chosen the parameter sets (J/hu,g) = (1, 1) and (1, \/2) largely because these 
points are commonly transected by the ground state curves of different approaches in the 
literature and thus permit the widest possible comparisons to be made. We have also chosen 
these parameters values because they are moderate, falling in a regime where no particular 
limiting theory is favored. Indeed, several results included in our tabulations (indicated by 
dagger symbols) are well-known asymptotic results not strictly valid in these intermediate 
scenarios. 

In making the selection of data to display in Tables [l| and [|, we have not attempted to be 
exhaustive, but to present a sampling of results from differing methodologies, highlighting 
more recent examples. With the exception of the DMFT values, all tabulated data have been 
computed directly by us, obtained from cited analytic formulae, or obtained directly from 
the original authors by private communication. The values shown for dynamical mean field 
theory have been measured from published figures, and are believed by us to fairly represent 
the published data to the number of significant digits displayed. 

Since not all of the approaches sampled in Tables [l] and |2| are variational calculations, 
we should note that a comparison of a variational energy E var with a lower non-variational 
energy E non (E var > E non ) in itself yields no conclusion regarding which value is "better" 
unless there exists independent proof that the non-variational energy lies above the ground 
state. In the opposite circumstance (E non > E var ), however, the variational principle alone 
assures that the variational energy is the better estimate. 
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Table 1: Ground State Energies as computed by various methods for (J/hu,g) = (1,1), 
(7, A) = (|, 1). The table is broken at the most likely location of the true ground state 
energy Values indicated by dagger symbols are obtained from limiting formulae that at 
these parameter values are clearly beyond their limits of validity 



Value Kuj] 


Type of Method 


Source 


- 1.73575 f 


AT 1 * 1 j_ * 11 1 

JNon-adiabatic small polaron 


TT 1 j M TH R 

Holstem,H Eq. |9j 


- 1.97 


QMC iV = 32/3 = lm = 10 


De Raedt and LagendijkErli-3 


nnnnn 4. 
- Z.UUUUU | 


Adiabatic strong coupling 


ljOgoIin,EJ t,q. |i(J| 




"»% 1 1 Wl IP d CC1 PQ "\ TCI VI QT1PT1 


T^"" ClInOQ l^Q Q A 11 V^T 1 ! T QTirl Toi T*/~\T1 ^ 

IVcLlUbcLivcLo, xiUUij, dllU. J. bllUIllbl — II — 1 


9 ^ 


OA/TP N — 39 /9 — ^ m — 39 

V^IVIO iv — OZ ^ — O ill — oz 


1 ^ O T-? O O/"] T" O TT /"I T O ITQT1 *^ 1 1 l^Pll 10 

l^t: JrtcLeClL all (J. J^cx^cI1Q1JKIjL_J 


- 2.40 


Dynamical Mean Field 


Ciuchi, et alB 


- 2.44721 


2 nd Qrder W CPT 




- 2.45611 


Merrifield variation 


Zhao, Brown, and Lindenberg^l 


- 2.46 


QMC N = 32 (3 = 20 m = 256 


De Raedt0 


- 2.46869 


Toyozawa variation 


Zhao, Brown, and LindenbergHI 


- 2.46931 


Global-Local variation 


Figure [1], this paper ^ ^ 


- 2.46968 


DMRG V = 32 


Jeckelmann and White0ll 


- 2.471 


Cluster diagonalization, N = 6 


Alexandrov, et a 1.0 BE 


- 3.08959 f 


2 nd order SCPT 


Marsiglio,il Stephan,@ Eq. UJ 



The most probable value of the true ground state energy can be sorted out from the data 
in Tables p] and |2| as follows: In view of the fact that the uncertainties in the best QMC values 
tabulated are of the order of one percent or greater ,0 the QMC values cannot be viewed as 
any more definitive than the Global-Local, DMRG, and cluster values that are all determined 
to higher precision and fall within the QMC uncertainties. Further, since the density matrix 
renormalization group method has been shown to have a variational realization,!!^ it would 
appear that the true ground state energy should lie below the DMRG values. On the other 
hand, as is discussed in greater detail in Section 4, cluster energies appear to converge upward 
with increasing cluster size, suggesting that the true bulk ground state energy lies above the 
N = 6 cluster values. Thus, each of the tables shown has been broken at the point where 
the best apparent upper bound to the ground state energy meets the best apparent lower 
bound. 

Turning our attention to the broader landscape of the system parameter space, Figure [l] 
shows the dependence of the global ground state energy E(0) on the coupling strength for 
assorted values of J/%uo from 1/4 to 5. The overall behavior of the ground state energy is 
to trend between two asymptotic g 2 dependences with differing coefficients and offsets. The 
crossover between these asymptotic trends occurs at the self trapping transition (g « gsr), 
which through analyses presented elsewhereM has been found to be accurately located by 
the relation g S r = 1 + \J J /Uoj. 

Using weak-coupling perturbation theory, one can show that the leading dependence of 



5 



Table 2: Ground State Energies as computed by various methods for (J/hu,g) = (1, y/2), 
(7, A) = (|,1). The table is broken at the most likely location of the true ground state 
energy. Values indicated by dagger symbols are obtained from limiting formulae that at 
these parameter values are clearly beyond their limits of validity 



T t 1 r-*- "1 

Value [nu>\ 


m pur 1 1 l 

type of Method 


Source 


- 2.27067 f 


Non-adiabatic small polaron 


HolsteinE Eq. | 


- 2.50000 f 


Adiabatic strong coupling 


Gogolm,EI Eq. 


- 2.51815 


bemiclassical variation 


Kalosakas, Aubry, and lsiromsl-ju 


- 2.73 


QMC N = 32p = lm = l0 


De Kaedt and LagendykEtlia 


- 2.8b 


QMC i\=62p = om = S2 


De Kaedt and LagendijkErli^ 


- 2.89 


Dynamical Mean Field 


Ciuchi, et al.@ 


- 2.89442 


2 nd Qrder W CPT 


Eq. | 


- 2.93301 


Merrifield variation 


this paper 


- 2.99172 


Toyozawa variation 


this paper 


- 2.99802 


Global-Local variation 


Figure [1], this paper 


- 2.99883 


DMRG iV = 32 


Jeckelmann and White0El 


- 3.000 


Cluster diagonalization, N = 6 


Alexandroy, et al.000 


- 3.00 


QMC N = 32 = 20 m = 256 


De Raedt0 


- 3.05279 f 


2 nd order SCPT 


Marsiglio,il Stephan,@ Eq. 



the ground state energy on the coupling constant is given for any J by 

g 2 Tiuj 



E(0) w -2 J 



1 + 4 J/hu 



This is, in fact, what we find to within numerical precision within the weak coupling regime 
of the Global-Local method (see Figure [I]). 

The weak-coupling estimate (El) is superior to the non-adiabatic small polaron esti- 
mate§§'@ 

E(0) « -g 2 huo - 2Je~ 92 (9) 

for most J provided the coupling strength is not too large (g < gsr), though the two 
approximations agree when both J/Tiu and g are small. When extrapolated into the adiabatic 
regime, however, @ is not even qualitatively sensible at weak coupling for J/hu > 1/2, and 
does not even approach the weak-coupling estimate until g ~ gsr- At strong coupling 
(9 ^ 9st)i (9) approaches —g 2 hu> for any J, which generally differs significantly from the 
true value of the ground state energy. Thus, it appears that there is nowhere in the adiabatic 
regime where the non-adiabatic small polaron approximation provides a meaningful estimate 
of the ground state energy. 

On the other hand, the adiabatic strong-coupling perturbation result (in one dimension)El 

£ (0)«-2JA(l + i L)=- 9 ^--£L (10) 
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Figure 1: The global ground state energy E(0) vs. g 2 for J/hou = 1/4, 1/2, 1, 2, 3, 4, 
and 5 in order from top to bottom (note that \im g ^ E(0) = —2 J). The straight dotted 
lines asymptoting each Global-Local curve at small g 2 are given by the second-order WCPT 
formula (||). The arched dashed lines asymptoting each Global-Local curve at large g 2 are 
given by the second-order SCPT formula (|10D. 

correctly describes the asymptotic behavior at strong coupling for g ^> g$T- Both this 
adiabatic result and the non-adiabatic result (§) can be recovered from the second-order 
strong-coupling perturbation result"' eI 

E(0) « -g 2 huo -2Je- g2 

-^J 2 e^ 2 [f(2g 2 )+f(g 2 )), (11) 

f(x) = Ei(x) - 7 - ln(x) , (12) 

in which Ei(x) is the exponential integral and 7 is the Euler constant. Examples of this 
strong-coupling result are plotted in Figure together with comparable Global-Local and 
weak-coupling results. The second order SCPT result is good for all g provided J/hu < 1/4, 
and for all J provided g 3> gsr, but breaks down rather dramatically otherwise. 

Thus, while both weak- and strong-coupling perturbation theories are clearly limited in 
scope, both are in quantitative agreement with results of the Global-Local method in the 
appropriate regimes, suggesting that a quite complete picture is available in the results of 
the Global-Local method augmented by the leading-orders of perturbation theory. Indeed, 
this blended approach proves to be broadly useful, as is borne out in results presented 
elsewhere r 
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4 Energy Bands 

The archetypical result of polaron theory in the strong coupling limit is the band form!' 

E(k) = {E- 2J) -E b - -B[l - cos(k)] . (13) 

The detailed dependence of the binding energy Eb and bandwidth B on J and g depend on 
regime, and it is through this dependence that changes in polaron structure are manifested in 
this limiting result. The binding energy and bandwidth are related quite simply in the non- 
adiabatic regime (E B ~ g 2 hu, B ~ 4Je -9 ), and less simply in the adiabatic regime. When 
this band form is valid, the polaron is heavily dressed by phonons and the polaron bandwidth 
B is small relative to both the bare electron bandwidth 4J and the bare phonon energy hu; 
i.e., B <C Min{4J, hu}. This narrowing of the polaron band and the related increase in the 
effective mass are commonly characterized as aspects of polaron band distortion. 

More generally, however, away from the strong coupling limit, polaron energy bands are 
not simple narrowed and shifted replicas of the bare band, but are more strongly distorted 
shapes whose non-sinusoidal dependence on k is a crucial reflection of polaron structure. For 
such bands, the polaron binding energy, effective mass, and polaron bandwidth no longer 
stand in any simple relationship to each other. 

In the limit g — > + of the adiabatic regime (J/hu > 1/4), one finds that the polaron 
energy band assumes a clipped form, 

E(k) = (E - 2 J) + 2 J[l - cos(k)] |«|<« c , (14) 
= (E — 2J)+hu \k\ > k c , (15) 

reflecting the difference in the character of polaron states above and below the wave vector 
k c (given by 2J[1 — cos(« c )] = hu) at which the bare electron energy band crosses into the 
one-phonon continuum. Although this crisply clipped band form is strictly valid only in the 
limit g — > + , its essential characteristics persist to non-trivial values of the electro n-phonon 
coupling strength. Figures |2] and [|, for example, demonstrate the persistence of the clipped 
band form and of k c as the relevant scale parameter for g of order unity. 

Figure |2] compares a complete Global-Local energv band with a comparable energy band 
as computed independently by the DMRG methocf3E§ for the case J/hu = 1, g = 1. As was 
the case in the ground state energy comparisons of Tables [IJ and 0, the DMRG energies are 
slightly lower than the GL energies, the difference ranging from about 0.015% at k = to 
about 0.46% at k = ir; no fitting or numerical adjustment of any kind was made to normalize 
these two results. 

The several surfaces completely detailing the electron-phonon structure underlying this 
energy band are shown in Figure |3|. 

These surfaces show polaron structure to be composed of multiple distinguishable and 
characteristic features and correlations. It is not our purpose here to analyze this internal 
structure in depth;B 001^0 however, a few outstanding features are worth noting. First, 
in all components of the solution, there exist distinguishable characters in the inner Brillouin 
zone (\k\ < k c ) and the outer Brillouin zone (\k\ > k c ). These distinct characters are least 
evident in the electron amplitudes a", and most evident in the primary phonon amplitude 
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Figure 2: Comparison of Global-Local (GL) energy band results with comparable results of 
the density matrix renormalization group (DMRG) method for J/fvuj = 1, g = 1. k c = 7r/3. 
DMRG data kindly provided by E. Jeckelmann.S. 

/3g. In the inner zone, phonon structure consists of phonon amplitudes correlated with the 
electronic component in a manner that is structured, but largely local in character. In the 
outer zone, phonon structure continues to have a localized component, but this component is 
dominated by a strong, momentum-rich, largely-delocalized component reflecting the strong 
influence of the one-phonon continuum on "clipped" polaron energy bands. Although these 
inner and outer zone features become muted with increasing coupling strength and become 
nearly indistinguishable above the self-trapping transition, their qualitative character is per- 
vasive. 

Figure ^ overlays the appropriate variational energy band computed by the Global-Local 
method upon data from several finite-cluster diagonalization results for the particular case 
J /hid = 1.25 and g = V0.5 • 1.25 ~ 0.79.... No fitting has been performed to achieve the 
impressive degree of agreement evident in these independent results. The variational results 
and the N = 20 cluster results agree to multiple significant digits across the entire Brillouin 
zone, the relative difference being approximately 0.01% at k = and 0.4% at k = tt. This 
degree of agreement is actually even better than it appears at first glance, as may be inferred 
from the convergence trends in the cluster data: 

Unlike variational bounds that converge downward toward the exact target energy, the 
trends evident in the data of Wellein and Fehske@'@ suggest that cluster energies converge 
upward with increasing cluster size. Figure |5] displays the ratio E^ =w (k) /\E^ 32 (k)\ for 
several k and N, where E^(k) is the band energy determined by Lanczos diagonalization 
on a cluster of N sites allowing up to M total phonons, and E^ 32 (k) is the Global-Local 
variational results computed for 32 sites as throughout this paper. The k = cluster energy 
is clearly well converged to a value slightly lower than, but nearly identical to the variational 
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Figure 3: Complete Global-Local solutions for the case J /Too = 1, g = 1 corresponding to 
the ground state energy tabulation in Table [I] and energy band presentation in Figure ^; a) 
Pq, b) Iq, c) Re{a£}, d) Im{afi}. k c = tt/3. 



energy. (That the k = cluster energy is nearly independent of iV suggests that the N = 6 
cluster energies0lll included in Tables |] and |^ are probably also well-converged.) Cluster 
energies are less well-converged away from the Brillouin zone center, where the evident trends 
suggest that the bulk limits of the fmite-K energies lie near the variational value. Thus, 
though in principle and through sufficiently comprehensive computation on sufficiently large 
clusters, cluster diagonalization methods should best our variational approach at any k, it 
appears that cluster sizes and retained phonon numbers must be significantly larger than 
those attempted to date in order improve significantly upon the Global-Local results for the 
polaron energy band as a whole. 

One conclusion that can be drawn from these comparisons (considering, for example, a 
fixed N) is that the outer Brillouin zone (|«| > k c ) displays a greater sensitivity to contribu- 
tions from higher phonon numbers than does the inner zone (\k\ < k c ). This is exactly what 
is to be expected from the K-dependence of the mean phonon number, which is typically 
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Figure 4: Comparison of our variational energy band with results of cluster diagonalizations 
for cluster sizes N = 10, 14, 18, and 20; J/Tiu = 1.25 and g = V0.5 ■ 1.25 = 0.79..., 
corresponding to (7, A) = (0.4,0.25). n c ~ 0.29527T. Cluster data kindly provided by G. 
WelleinH 

considerably smaller in the inner zone than in the outer zone (see, e.g., Figure 1 of Paper 
110). 

Owing to ultimate limitations of computing time and physical limitations of data storage, 
cluster diagonalizations generally are limited by the maximum dimension of the truncated 
Hilbert space that can be addressed by a particular computer. For a one-electron problem 
on a lattice of N sites, the dimension of the electronic subpace is D e i = N and the dimension 
D p h of a truncated phonon subspace containing at most M phonons is given by D p h = 
(N + M)\/N\M\. Calculations have been reported for various combinations of N and M from 
N = 2 to N = 20 and M up to 50, each such calculation striking a unique balance between N 
and M consistent with computational resources. Table ||] presents values conveying the scale 
of the diagonalization problem for clusters sizes and phonon numbers up to N = M = 32. 
Clearly, the scenarios represented by the lower right half of the tabulated values lie beyond 
the scope of machine diagonalization. 

Implementation of the density matrix renormalization group method also involves a trun- 
cation of the total quantum Hilbert space; however, the limitations imposed by such trun- 
cations appear to be less onerous than in the case of cluster diagonalization, permitting 
relatively large calculations to be accomplished in reasonable time on desktop platforms, 
much as the global local method. 

By contrast, the "size" of a Global-Local variation does not scale directly with phonon 
number, since the nature of the Global-Local approximation is not to truncate the phonon 
Hilbert space but to characterize the phonon distribution for M = 00 through a relatively 
small number of parameters. For example, in all the calculations presented in this paper, 
only 3N independent variational parameters were required to flexibly represent the polaron 
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Figure 5: Cluster energies for k = 0, 7r/2, and 7r and cluster sizes iV = 10 to 20 for 
J/Tiu = 1.25 and g = V0.5 • 1.25 ~ 0.79.... Each cluster energy is divided by the absolute 
value of our variational energy appropriate to that particular k, so that —1 indicates identical 
results under the two methods. Cluster data kindly provided by G. Wellein.il 

state to high precision at any k (96 for |«| 6 (0, it), 48 for |«| =0, 7r) with no restriction on 
phonon numbers, vastly fewer than comparable cluster diagonalizations. 

It is also the case that the Global-Local method is not at its best for small clusters, 
but improves in quality with increasing N (though not without limit). This improvement is 
realized because the total phonon state is represented by a superposition of N distinct phonon 
coherent states; such a superposition becomes more flexible as the number of coherent states 
in the superposition increases. 

Our own comparisons with 6-cluster ground state energies in Tables [1] and proved 
mutually favorable, as has our comparison with the N = 10 — 20 cluster ground states as 
computed by Wellein and Fehske.B'H However, we are more generally concerned with the 
complete band structure of the polaron, which is not well captured by calculations on small 
clusters as can be seen from the following: 

From weak-coupling nearly up to the self-trapping transition, polaron bands exhibit dis- 
tinct character for wave vectors above and below k c that must be resolved if one is to 
faithfully represent polaron structure. This underscores the importance in any finite-lattice 
calculation of maintaining the lattice size N in relation to the tunneling matrix element J 
such that at least one finite-K point is sampled on both the high and low sides of k c . For 
a given N, this implies that only for J G {J mm-, J max) can both the inner and outer polaron 
band structure be characterized independently, where 

2[1 — cos(27r/iV)J 
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Table 3: Cluster Facts 



liln^t pr 

LLO uv^l 


7 • Ifiuj 


7 Ifiuj 


D< + M — 1 


D+ 4 M — 1 6 

± -'tot i lvl — xu 


D + * M = 32 

^toti lvl — J* 


2 sites 


no vaiue 


no vaiue 




ouo 


1 1 99 


4 sites 


0.5000 


0.500 


4.0 x 10 3 


1.9 x 10 4 


2.4 x 10 5 


6 sites 


0.3333 


1.000 


4.8 x 10 4 


4.5 x 10 5 


1.7 x 10 7 


10 sites 


0.2764 


2.618 


1.8 x 10 6 


5.3 x 10 7 


1.5 x 10 10 


16 sites 


0.2599 


6.569 


8.5 x 10 7 


9.6 x 10 9 


3.6 x 10 13 


20 sites 


0.2563 


10.22 


6.0 x 10 s 


1.5 x 10 11 


2.5 x 10 15 


32 sites 


0.2524 


26.02 


4.7 x 10 10 


7.2 x 10 13 


5.9 x 10 19 



mm - 2[l-cos((iV-2)7r/JV)] ' 1 J 

Most crucially, for J > J maa; , it is not possible to resolve the parabolic band bottom that 
is characteristic of small k from which to extract a meaningful effective mass. To illustrate 
the gravity of this restriction, the values of J m i n and J max so defined have been presented in 
Table for each of the cluster sizes there considered. It is clear from these values that cluster 
size imposes a significant limitation on one's ability to describe polaron band structure over 
a meaningful interval of J. 



5 Effective Mass 

The ground state energy as considered in Section 3 permitted a fair comparison of many 
different approaches in part because in many cases approximate methodologies are at their 
"best" when applied to the global ground state; excited states almost universally pose greater 
challenges. In energy band theory, the minimal excursion beyond the ground state is con- 
tained in the effective mass, since in principle the determination of the effective mass requires 
knowledge of only an infinitesimal excitation to finite k. In this section we compare effective 
masses as yielded by a number of different approaches. 

Our effective mass computations were based on the formula 

m* _ 2J 

m ° -ay 

using a discrete representation of the k derivative at the Brillouin zone center. We note that 
the m used in our calculations was obtained by computing its value using the same discrete 
differentiation as was applied to compute m* rather than using the limiting iV — > oo value. 
Not only does this minimize any dependence of the computed effective mass ratio on lattice 
size, but is technically necessary in order to properly normalize our results; e.g., to exactly 
recover the limit lim 9 ^ m*/m = 1. For all cases presented, J < J ma x according to ( [TTD for 
lattices of 32 sites, so that our discrete differentiation yields the physically-meaningful value. 

In Figure |^, we first consider the simultaneous comparison of our Global-Local effective 
masses ( actually, ln(m*/mo) ) with comparable DMRG results and with second-order strong- 
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coupling perturbation theory as contained in the relation0@ 

m* e g2 

^ = 1 + (4J/Me- ff W) ' 

for which f(x) has been defined in fll2|) . 




Figure 6: Effective mass ratio. Bold curves: Global-Local results for J/Tuo = 1/4, 1, and 5. 
Faint curves: Second-order strong-coupling perturbation theory (see Eq. [19]). Points: Den- 
sity matrix renormalization group (DMRG) results as kindly provided by E. Jeckelmannilil. 

The Global-Local and DMRG results agree quite well, especially considering that these 
particular DMRG data do not represent actual finite-K calculations, but estimations based 
on DMRG ground state data.Sii Given the very favorable energy band comparison in 
Figure ||, it is reasonable to expect that a future comparison based on direct excited-state 
DMRG calculations would show even better agreement (e.g., see Figure |2|). What deviations 
there are between the GL and DMRG results are greatest above the self-trapping transition 
(g > gsr)] it is unclear at present whether this is a discrepancy of lasting significance. 

On the other hand, comparison of both GL and DMRG effective masses with the results 
of second-order SCPT are far less favorable. Agreement is excellent for J /Two = 1/4 and 
g > gsr'-, however, deviations appear even at this small J value for g < g$T- At J/Tioj = 1, 
it is evident that both the GL and DMRG masses asymptote to the SCPT result, but that 
this convergence of results does not materialize until well above the self-trapping transition 
(g ^ 9st)- At J/TujJ = 5, the disagreement between the SCPT result and both the GL 
and DMRG masses is so severe as to render the SCPT estimation useless. We note that 
although we have explicitly compared only selected J values, it can be safely inferred that 
the second-order SCPT mass estimation ceases to be relevant in practical terms by the time 
J/fiuj ~ 2.0 



14 



In Figure ^ we make a more narrow comparison with a broader range of approaches. 
Here, we again compare GL and DMRG results, this time including as wellS 

i) data from direct quantum Monte Carlo calculations of the polaron mass,00 

ii) the second-order SCPT estimate as contained in (191), 
hi) the second-order WCPT estimate as contained inE3 

mo _ g 2 (l + 2J/hu) 

m* (1 + 4J/M 3/2 



and iv) the weak-coupling Migdal estimate as contained in0.ll 



m* 2 (l + 2J/huj) 



i+ ^r,::^ . (2i) 



It is clear that all the presented results except the SCPT are mutually consistent at suf- 
ficiently weak coupling; however, divergences between results appear quickly with increasing 
coupling strength, such that only the GL, DMRG, and WCPT masses remain mutually 
consistent to nontrivial coupling. 

We note that the agreement between the GL, DMRG, and WCPT masses is actually 
better than it may appear. As in Figure || the DMRG data indicated by diamond symbols 
is the result of an estimation procedure based on the global DMRG ground state only. The 
near-perfect agreement between the GL and DMRG energy bands as shown in Figure ||] 
assures that the effective masses computed directly from the complete energy bands are 
essentially identical. The same high degree of agreement can be inferred from the comparison 
of GL and cluster diagonalization bands in Figure [|. The small apparent discrepancy between 
the indirect DMRG estimate and the GL and SCPT results is thus essentially eliminated, 
showing GL, DMRG, and WCPT to be in essentially complete agreement up to at least 
9 ~ I- 

At strong coupling, the GL, DMRG, and SCPT masses are again mutually consistent, 
while each of the remaining results presented deviates significantly. 



6 Conclusion 

In this paper, we have presented results for the Holstein molecular crystal model in one 
space dimension as determined by the Global-Local variational method, including complete 
polaron energy bands, ground state energies, and effective masses. We have juxtaposed our 
results with specific comparable results of numerous other methodologies of current interest, 
including quantum Monte Carlo, cluster diagonalization, dynamical mean field theory, den- 
sity matrix renormalization group, semiclassical analysis, weak-coupling perturbation theory, 
and strong-coupling perturbation theory. 

Through these comparisons, we have been able to conclude that 

i) perturbation theory, while definitive in limits, is clearly superceded by more accurate 
methods in the intermediate-coupling regime, 

ii) semiclassical analysis, while consistent with perturbation theory in the strong-coupling 
limit, does not improve significantly over perturbation theory away from this limit, 
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Figure 7: Inverse effective mass ratio for J/hu = 1. Solid curve: Global-Local result. 
Diamond symbols: DMRG results as kindly provided by E. JeckelmannE3c2l. Bullet symbols: 
QMC results as kindly provided by P. Kornilovit Dotted curve: Weak-coupling Migdal 

approximation (see Eq. |2T| ). Dashed curve: Second-order strong-coupling perturbation 
theory (see Eq. [19]). Dot-dashed curve: Second-order weak-coupling perturbation theory 
(see Eq. E3). 



iii) quantum Monte Carlo for the ground state energy is consistent with its best com- 
petitors, but owing to limitations of precision is in the present context more confirmatory of 
other approaches than determinative in itself, 

iv) quantum Monte Carlo for the effective mass is consistent with more accurate ap- 
proaches at weak coupling, but at present is inconsistent with the best estimates at strong 
coupling, 

v) dynamical mean field theory, while exceptionally complete qualitatively, does not 
deliver superior quantitative accuracy in the regimes considered in this paper. 

On the positive side, very favorable comparisons were found between the Global-Local 
method, the density matrix renormalization group method, and cluster diagonalization meth- 
ods. Of these, however, only the GL and DMRG methods could be compared in a broad way. 
Both GL and DMRG were found to be consistent with perturbation theory at both weak 
and strong coupling, and in the weak-coupling case, at least, consistent to higher coupling 
strengths than any other method. Where GL and DMRG depart from perturbation theoretic 
results, in the intermediate coupling regime, GL and DMRG remain mutually consistent to 
an impressive degree, and where direct comparison has been possible, consistent with cluster 
diagonalization. 

We take these results as a whole to confirm the Global-Local variational method as being 
highly accurate over a wide range of the polaron parameter space, from the non-adiabatic 
limit to the extremes of high adiabaticity, from weak coupling through intermediate coupling 
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to strong coupling. 

In succeeding works,i~i we shall present a comprehensive analysis of polaron structure 
and properties as determined in the Global-Local method over essentially the whole of the 
polaron parameter space. 
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Value [Huj] 


Type of Method 


Source 


- 1.73575 f 


Non-adiabatic small polaron 


HolsteinJ Eq. g 


- 1.97 


QMC N = 32 p = 1 m = 10 


De Raedt and Lagendijk!0 


- 2.00000 f 


Adiabatic strong coupling 


Gogolin,i Eq. ^ 


- 2.08614 


Semiclassical variation 


Kalosakas, Aubry, and TsironisS'S 


- 2.35 


QMC iV = 32/3 = 5m = 32 


De Raedt and Lagendijk§0 


- 2.40 


Dynamical Mean Field 


Ciuchi, et alB 


- 2.44721 


2 nd Qrder W CPT 




- 2.45611 


Merrifield variation 


Zhao, Brown, and Lindenbergti3 


- 2.46 


QMC N = 32 p = 20 m = 256 


De RaedtS 


- 2.46869 


Toyozawa variation 


Zhao, Brown, and LindenbergH3 


- 2.46931 


Global-Local variation 


Figure [l], this paper ^ ^ 


- 2.46968 


DMRG N = 32 


Jeckelmann and White0il 


- 2.471 


Cluster diagonalization, N = 6 


Alexandre^ et al.0@0 


- 3.08959 f 


2 nd order SCPT 


MarsiglioP StephanJH Eq. y 
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Value [huj] 


Type of Method 


Source 


- 2.27067 f 


Non-adiabatic small polaron 


Holstein^ Eq. 




- 2.50000 f 


Adiabatic strong coupling 


Gogolin,y Eq. 


m 


- 2.51815 


Semiclassical variation 


Kalosakas, Aubry, and TsironisBS 


- 2.73 


QMC N = 32 (3 = 1 m = 10 


De Raedt and Lagendijk§0 


- 2.86 


QMC iV = 32/3 = 5m = 32 


De Raedt and Lagendijki0 


- 2.89 


Dynamical Mean Field 


Ciuchi, et alB 


- 2.89442 


2 nd Qrder W CPT 


Eq. | 


- 2.93301 


Merrifield variation 


this paper 


- 2.99172 


Toyozawa variation 


this paper 


- 2.99802 


Global-Local variation 


Figure [l], this paper 


- 2.99883 


DMRG N = 32 


Jeckelmann and Whitellil 


- 3.000 


Cluster diagonalization, N = 6 


Alexandrov, et al0EE0 


- 3.00 


QMC iV = 32 p = 20 m = 256 


De RaedtS 


- 3.05279 f 


2 nd Qrder g C p T 


Marsiglio,il Stephan,@ Eq. y 
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Cluster 






D tot , M = 10 


D tot ,M = 16 


D tot) M = 32 


2 sites 


no value 


no value 


132 


306 


1122 


4 sites 


0.5000 


0.500 


4.0 x 10 3 


1.9 x 10 4 


2.4 x 10 5 


6 sites 


0.3333 


1.000 


4.8 x 10 4 


4.5 x 10 5 


1.7 x 10 7 


10 sites 


0.2764 


2.618 


1.8 x 10 6 


5.3 x 10 7 


1.5 x 10 1U 


16 sites 


0.2599 


6.569 


8.5 x 10 7 


9.6 x 10 9 


3.6 x 10 13 


20 sites 


0.2563 


10.22 


6.0 x 10 8 


1.5 x 10 11 


2.5 x 10 15 


32 sites 


0.2524 


26.02 


4.7 x 10 1U 


7.2 x 10 13 


5.9 x 10 19 
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